In [34]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['gamma', 'test']
`%matplotlib` prevents importing * from pylab and numpy

Traffic Models with Discontinuous and Non-convex Flux


Qi Guo, Peng Zheng

University of Washington

Overview

  • Freeway driving
    • Introduction
    • Method of mollification
    • Exact solution, CLAWPACK solution and car-following model
  • Night-time driving
    • Introduction
    • Clustering beharvior

LWR Model

In 1955, Lighthill, Whitham and Richards proposed the a model, now known as LWR model for traffic flow problem. This model is a macroscopic traffic flow model.

$$\rho_t+f(\rho)_x=0,$$

where $$f(\rho)=\rho v(\rho).$$

Freeway Driving Model

From the empirical data collected from freeway, the velocity can be modeled as:

$$v(\rho)=\left\{\begin{array}{lll} 1,&\text{if }0\leq\rho<\rho_m,\\ \gamma(\frac{1}{\rho}-1),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$

Hence

$$f(\rho)=\left\{\begin{array}{lll} \rho,&\text{if }0\leq\rho<\rho_m,\\ \gamma(1-\rho),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$

For this problem we use $\rho_m=0.5$ and $\gamma=0.5$.

In [35]:
rhom = 0.5
gamma = 0.5
def flux_freeway_origin(rho):
    n = rho.size
    f = zeros(n)
    for i in range(n):
        if rho[i] >= 0. and rho[i] < rhom:
            f[i] = rho[i]
        elif rho[i] >= rhom and rho[i] <= 1.:
            f[i] = gamma*(1.-rho[i])
        else:
            print "Invalid density value!"
    return f
In [36]:
rho1 = linspace(0.,rhom-0.01,50)
rho2 = linspace(rhom,1.,51)
plot(rho1,flux_freeway_origin(rho1),'b')
plot(rho2,flux_freeway_origin(rho2),'b')
plot([rhom,rhom],[rhom,gamma*(1.-rhom)],'--b')
ylim(0.,0.5)
xlabel('rho')
ylabel('$f$')
title('Original flux fucntion')
Out[36]:
<matplotlib.text.Text at 0x1089f9ad0>

Method of Mollification

The discontinuity can result the invalidation of the technique we learned from the class, since the mehod involves the deriavtive of the flux function. One method to redeem this is to use mollifier smooth the flux function. We use the mollifier as follows:

$$\eta_\epsilon(s)=\frac{1}{\epsilon}\eta\left(\frac{s}{\epsilon}\right),$$

where

$$\eta(x)=\left\{\begin{array}{lll} C\exp\left(\frac{1}{x^2-1}\right)&\text{if }-1< x < 1,\\ 0&\text{otherwise},\end{array}\right.$$

Here $C$ is a constant, $C\approx2.2522836.$

The modified flux function is the convolution of the mollifier and original function:

$$f_\epsilon=\eta_\epsilon\ast f,$$

so

$$f_\epsilon(\rho)=\rho+(\gamma(1-\rho)-\rho)\int_{\rho_m-\epsilon}^\rho\eta_\epsilon(s-\rho_m)ds,$$

with $0<\epsilon\ll1$.

It has been proven that, if $\eta_\epsilon(x)\rightarrow\delta(x)$ as $\epsilon\rightarrow 0$, then $f_\epsilon\rightarrow f$.

In the project, we select $\epsilon=0.01$.

In [37]:
def eta(x,eps):
    C = 2.25228362104358101049978125556
    eta = where(abs(x)<eps,C*exp(eps**2/(x**2-eps**2))/eps,0.)
    return eta
def flux_freeway_modified(rho,eps):
    n = rho.size
    f = zeros(n)
    
    for i in range(n):
        N = 500*(int((rho[i]-rhom+eps)/eps) + 1) # mesh cell for integration
        dx = (rho[i]-rhom+eps)/N # mesh size for integration
        intl = 0.
        for j in range(N):
            intl = intl + dx*eta((j+0.5)*dx-eps,eps)
        f[i] = rho[i] + (gamma*(1.-rho[i])-rho[i])*intl
    return f
In [38]:
rho = linspace(0.,1.,101)
hold
plot(rho,flux_freeway_modified(rho,0.01),'m',label='$\epsilon=0.01$')
plot(rho,flux_freeway_modified(rho,0.1),'b',label='$\epsilon=0.1$')
ylim(0.,0.5)
xlabel('$rho$')
ylabel('$f_\epsilon$')
legend()
/Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:11: RuntimeWarning: divide by zero encountered in double_scalars
/Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: divide by zero encountered in double_scalars
  app.launch_new_instance()
/Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: overflow encountered in exp
  app.launch_new_instance()
Out[38]:
<matplotlib.legend.Legend at 0x107a68c50>

Exact Solution, Numerical Solution and Car-Following Model

Overview

Explicit Method

To determine the weak solution to a non-convex scalar conservation law, in this problem, we will use Olenik entropy condition:

$$ \frac{f(\rho)-f(\rho_l)}{\rho-\rho_l}\geq \frac{f(\rho_l)-f(\rho_r)}{\rho_l-\rho_r}\geq \frac{f(\rho)-f(\rho_r)}{\rho-\rho_r} $$

for all $\rho$ between $\rho_l$ and $\rho_r$.

The entropy-satisfying solution to this problem can be determined form the graph of $f_\epsilon(\rho)$. This method is called convex hull.

Numerical Method

We use the 2nd order Lax-Wendroff method with Van Leer limiter.

In [39]:
import glob
from matplotlib import image
from clawpack.visclaw.JSAnimation import IPython_display
from matplotlib import animation

def init():
    im.set_data(image.imread(filenames[0]))
    return im,

def animate(i):
    image_i=image.imread(filenames[i])
    im.set_data(image_i)
    return im,

Car-Following Model

We only focus on each car $x_n(t)$.

The density between two car is the reciprocal of the distance: $$\rho_{n+1}(t)=\frac{1}{x_{n}(t)-x_{n+1}(t)};$$

The velocity of each car only depends on the distance: $$v_{n}(t)=v(\rho_{n}(t));$$

The position is updated by: $$x_{n}(t+\Delta t)=x_{n}(t)+v_{n}(t)\Delta t.$$

Based on these assumptions, we can build up our own ODE Solver.

In [40]:
# ODE Solver Initialization

def initialization(N,rhol,rhor):
    return linspace(0.,-(N-1.)/rhol,N)
In [41]:
# Main Function of ODE Solver

def ode_solver_freeway(N,tfinal,output_number,rhol,rhor):
    from copy import copy
    Dt = tfinal/output_number
    dt = min(0.01,Dt)
    x = zeros((output_number+1,N))
    rho = copy(x)
    x[0] = initialization(N,rhol,rhor)
    rho[0][1:N+1] = -1./diff(x[0])
    rho[0][0] = rhor
    m = int(Dt/dt)
    xs ,rhos = copy(x[0]), copy(rho[0])
    for i in range(1,output_number+1):
        for j in range(m):
            v = flux_freeway_origin(rhos)/rhos
            xs = xs + v*dt
            rhos[1:N+1] = -1./diff(xs)
            rhos[0] = rhor
        v = flux_freeway_origin(rhos)/rhos
        xs = xs + v*(Dt-dt*m)
        rhos[1:N+1] = -1./diff(xs)
        rhos[0] = rhor
        x[i], rho[i] = copy(xs), copy(rhos)
    return x, rho
In [110]:
N = 80
tfinal = 40.
output_number = 80

Case 1: $\rho_r<\rho_m<\rho_l$

Exact Solution

We construct the smallest convex hull of the set $\{(\rho,y):\rho_r<\rho<\rho_l,\ \text{and}\ y\leq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:

$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq t,\\ \rho_r \ \ \ \text{if}\ \ x > t, \end{array} \right.$$

where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$, since the shock and characteristic speeds are equal.

In the numerical simulation, the initial condition we are using is:

$$\rho_l=0.9, \rho_r=0.2.$$

PDE Solver

In [43]:
%%system

make data ql=0.9 qr=0.2
make .plots
Out[43]:
["make: `.data' is up to date.", "make: `.plots' is up to date."]
In [44]:
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[44]:


Once Loop Reflect

ODE Solver

In [112]:
rhol = 0.9
rhor = 0.2
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
In [46]:
%%system
mkdir _plots_ode_freeway1
Out[46]:
['mkdir: _plots_ode_freeway1: File exists']
In [113]:
Dt = tfinal/output_number
for i in range(output_number+1):
    clf()
    figure(figsize=(12,5))
    subplot(1,2,1)
    for j in range(0,N,4):
        plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
    plot(x[i][0:N:4],Dt*i*ones(N/4),'ob')
    plot([-N/rhol-10,tfinal+10],[Dt*i,Dt*i],'--r');
    axis([-N/rhol-10,tfinal+10,0,tfinal+1])
    title('Trajectory of Cars')
    xlabel('$x$')
    ylabel('$t$')
    
    subplot(1,2,2)
    plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
    axis([-N/rhol-10,tfinal+10,-0.1,1.1])
    title('Density')
    xlabel('$x$')
    ylabel('rho')
    
    fname = '_plots_ode_freeway1/frame%04i'%i + 'fig1' + '.png'
    savefig(fname)
<matplotlib.figure.Figure at 0x111748050>
<matplotlib.figure.Figure at 0x11179e0d0>
<matplotlib.figure.Figure at 0x10d71bd90>
<matplotlib.figure.Figure at 0x11a7876d0>
<matplotlib.figure.Figure at 0x124d96d90>
<matplotlib.figure.Figure at 0x128d09f50>
<matplotlib.figure.Figure at 0x1243c7610>
<matplotlib.figure.Figure at 0x1036b6f10>
<matplotlib.figure.Figure at 0x11166f650>
<matplotlib.figure.Figure at 0x10fba60d0>
<matplotlib.figure.Figure at 0x10e121dd0>
<matplotlib.figure.Figure at 0x1140dcd50>
<matplotlib.figure.Figure at 0x11c3c17d0>
<matplotlib.figure.Figure at 0x113f38790>
<matplotlib.figure.Figure at 0x113f02490>
<matplotlib.figure.Figure at 0x11a73f4d0>
<matplotlib.figure.Figure at 0x10e4d6c50>
<matplotlib.figure.Figure at 0x10cc99a50>
<matplotlib.figure.Figure at 0x128cffc50>
<matplotlib.figure.Figure at 0x11ba0e9d0>
<matplotlib.figure.Figure at 0x1116ae810>
<matplotlib.figure.Figure at 0x10d6a5490>
<matplotlib.figure.Figure at 0x11ba92510>
<matplotlib.figure.Figure at 0x11289fc50>
<matplotlib.figure.Figure at 0x10f63e7d0>
<matplotlib.figure.Figure at 0x119be1810>
<matplotlib.figure.Figure at 0x1160e1bd0>
<matplotlib.figure.Figure at 0x10ec68090>
<matplotlib.figure.Figure at 0x128d2e650>
<matplotlib.figure.Figure at 0x1185a57d0>
<matplotlib.figure.Figure at 0x11a6cb390>
<matplotlib.figure.Figure at 0x11ef96d50>
<matplotlib.figure.Figure at 0x11c61f990>
<matplotlib.figure.Figure at 0x113f386d0>
<matplotlib.figure.Figure at 0x113e4fe90>
<matplotlib.figure.Figure at 0x1185d6cd0>
<matplotlib.figure.Figure at 0x11a77aad0>
<matplotlib.figure.Figure at 0x11c6b4b90>
<matplotlib.figure.Figure at 0x1243cb210>
<matplotlib.figure.Figure at 0x124d96d50>
<matplotlib.figure.Figure at 0x120b04c90>
<matplotlib.figure.Figure at 0x10e0f29d0>
<matplotlib.figure.Figure at 0x10ff8d190>
<matplotlib.figure.Figure at 0x10e11e890>
<matplotlib.figure.Figure at 0x11e268cd0>
<matplotlib.figure.Figure at 0x113e53450>
<matplotlib.figure.Figure at 0x10ce11790>
<matplotlib.figure.Figure at 0x113f0d190>
<matplotlib.figure.Figure at 0x11ef7a390>
<matplotlib.figure.Figure at 0x10d7cf6d0>
<matplotlib.figure.Figure at 0x113e53090>
<matplotlib.figure.Figure at 0x11ba0e450>
<matplotlib.figure.Figure at 0x10ec40710>
<matplotlib.figure.Figure at 0x11a6cbf10>
<matplotlib.figure.Figure at 0x11a6b6290>
<matplotlib.figure.Figure at 0x119c15290>
<matplotlib.figure.Figure at 0x1128d8110>
<matplotlib.figure.Figure at 0x128ee5410>
<matplotlib.figure.Figure at 0x1128a32d0>
<matplotlib.figure.Figure at 0x11c3c1e10>
<matplotlib.figure.Figure at 0x10eca8090>
<matplotlib.figure.Figure at 0x10f61f610>
<matplotlib.figure.Figure at 0x10e4eb610>
<matplotlib.figure.Figure at 0x127572d10>
<matplotlib.figure.Figure at 0x1243d6190>
<matplotlib.figure.Figure at 0x1185a5350>
<matplotlib.figure.Figure at 0x124d7c510>
<matplotlib.figure.Figure at 0x10e4e1450>
<matplotlib.figure.Figure at 0x1243b0e90>
<matplotlib.figure.Figure at 0x114a3d510>
<matplotlib.figure.Figure at 0x127522190>
<matplotlib.figure.Figure at 0x111797fd0>
<matplotlib.figure.Figure at 0x10e13b350>
<matplotlib.figure.Figure at 0x10ff99c10>
<matplotlib.figure.Figure at 0x10f739e50>
<matplotlib.figure.Figure at 0x11c6a8050>
<matplotlib.figure.Figure at 0x11169bf50>
<matplotlib.figure.Figure at 0x113e40950>
<matplotlib.figure.Figure at 0x10d725f10>
<matplotlib.figure.Figure at 0x10ffa7190>
<matplotlib.figure.Figure at 0x10f6f1e50>
In [114]:
figno = 1
fname = '_plots_ode_freeway1/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(12, 8),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[114]:


Once Loop Reflect

Case 2: $\gamma/(1+\gamma)<\rho_l<\rho_m<\rho_r$

Exact Solution

We construct the smallest convex hull of the set $\{(\rho,y):\rho_l<\rho<\rho_r,\ \text{and}\ y\geq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:

$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq -\gamma t,\\ \rho_r \ \ \ \text{if}\ \ x > -\gamma t, \end{array} \right.$$

where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$.

In the numerical simulation, the initial condition we are using is:

$$\rho_l=0.4, \rho_r=0.9.$$

PDE Solver

In [54]:
%%system

make data ql=0.4 qr=0.9
make .plots
Out[54]:
['rm -f .data',
 'Wrote data to initial_data.txt',
 'python setrun.py              classic                  ',
 'touch .data',
 '/Library/Developer/CommandLineTools/usr/bin/make output',
 'rm -f .output',
 'python /Users/Peng/Documents/AMATH574/clawpack/clawutil/src/python/clawutil/runclaw.py xclaw                          _output                     \\',
 '\tTrue                    False                     . False False None',
 'Reading data file: claw.data   ',
 '         first 5 lines are comments and will be skipped',
 ' running...',
 '  ',
 'Reading data file: setprob.data',
 '         first 5 lines are comments and will be skipped',
 'CLAW1EZ: Frame    1 output files done at time t =  0.5000D-01',
 '',
 'CLAW1EZ: Frame    2 output files done at time t =  0.1000D+00',
 '',
 'CLAW1EZ: Frame    3 output files done at time t =  0.1500D+00',
 '',
 'CLAW1EZ: Frame    4 output files done at time t =  0.2000D+00',
 '',
 'CLAW1EZ: Frame    5 output files done at time t =  0.2500D+00',
 '',
 'CLAW1EZ: Frame    6 output files done at time t =  0.3000D+00',
 '',
 'CLAW1EZ: Frame    7 output files done at time t =  0.3500D+00',
 '',
 'CLAW1EZ: Frame    8 output files done at time t =  0.4000D+00',
 '',
 'CLAW1EZ: Frame    9 output files done at time t =  0.4500D+00',
 '',
 'CLAW1EZ: Frame   10 output files done at time t =  0.5000D+00',
 '',
 'CLAW1EZ: Frame   11 output files done at time t =  0.5500D+00',
 '',
 'CLAW1EZ: Frame   12 output files done at time t =  0.6000D+00',
 '',
 'CLAW1EZ: Frame   13 output files done at time t =  0.6500D+00',
 '',
 'CLAW1EZ: Frame   14 output files done at time t =  0.7000D+00',
 '',
 'CLAW1EZ: Frame   15 output files done at time t =  0.7500D+00',
 '',
 'CLAW1EZ: Frame   16 output files done at time t =  0.8000D+00',
 '',
 'CLAW1EZ: Frame   17 output files done at time t =  0.8500D+00',
 '',
 'CLAW1EZ: Frame   18 output files done at time t =  0.9000D+00',
 '',
 'CLAW1EZ: Frame   19 output files done at time t =  0.9500D+00',
 '',
 'CLAW1EZ: Frame   20 output files done at time t =  0.1000D+01',
 '',
 'CLAW1EZ: Frame   21 output files done at time t =  0.1050D+01',
 '',
 'CLAW1EZ: Frame   22 output files done at time t =  0.1100D+01',
 '',
 'CLAW1EZ: Frame   23 output files done at time t =  0.1150D+01',
 '',
 'CLAW1EZ: Frame   24 output files done at time t =  0.1200D+01',
 '',
 'CLAW1EZ: Frame   25 output files done at time t =  0.1250D+01',
 '',
 'CLAW1EZ: Frame   26 output files done at time t =  0.1300D+01',
 '',
 'CLAW1EZ: Frame   27 output files done at time t =  0.1350D+01',
 '',
 'CLAW1EZ: Frame   28 output files done at time t =  0.1400D+01',
 '',
 'CLAW1EZ: Frame   29 output files done at time t =  0.1450D+01',
 '',
 'CLAW1EZ: Frame   30 output files done at time t =  0.1500D+01',
 '',
 'CLAW1EZ: Frame   31 output files done at time t =  0.1550D+01',
 '',
 'CLAW1EZ: Frame   32 output files done at time t =  0.1600D+01',
 '',
 'CLAW1EZ: Frame   33 output files done at time t =  0.1650D+01',
 '',
 'CLAW1EZ: Frame   34 output files done at time t =  0.1700D+01',
 '',
 'CLAW1EZ: Frame   35 output files done at time t =  0.1750D+01',
 '',
 'CLAW1EZ: Frame   36 output files done at time t =  0.1800D+01',
 '',
 'CLAW1EZ: Frame   37 output files done at time t =  0.1850D+01',
 '',
 'CLAW1EZ: Frame   38 output files done at time t =  0.1900D+01',
 '',
 'CLAW1EZ: Frame   39 output files done at time t =  0.1950D+01',
 '',
 'CLAW1EZ: Frame   40 output files done at time t =  0.2000D+01',
 '',
 'Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL',
 '==> runclaw: Will take data from  /Users/Peng/Documents/AMATH574/am574-group07/presentation',
 '==> runclaw: Will write output to  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '==> runclaw: Removing all old fort files in  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '',
 '==> Running with command:',
 '    /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw',
 '',
 '==> runclaw: Finished executing',
 '',
 '==> runclaw: Done executing /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw via clawutil.runclaw.py',
 '==> runclaw: Output is in  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '/Library/Developer/CommandLineTools/usr/bin/make plots',
 'rm -f .plots',
 'python /Users/Peng/Documents/AMATH574/clawpack/visclaw/src/python/visclaw/plotclaw.py _output                     _plots                     setplot.py           ',
 'Importing setplot.setplot from /Users/Peng/Documents/AMATH574/am574-group07/presentation.',
 'Executed setplot successfully',
 'Will plot 41 frames numbered: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]',
 'Will make 1 figure(s) for each frame, numbered:  [1]',
 '',
 '-----------------------------------',
 '',
 '',
 'Creating html pages for figures...',
 '',
 "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ",
 '    already exists, files may be overwritten ',
 'Now making png files for all figures...',
 '    Reading  Frame 0 at t = 0  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 0 at time t = 0.0',
 '    Reading  Frame 1 at t = 0.05  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 1 at time t = 0.05',
 '    Reading  Frame 2 at t = 0.1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 2 at time t = 0.1',
 '    Reading  Frame 3 at t = 0.15  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 3 at time t = 0.15',
 '    Reading  Frame 4 at t = 0.2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 4 at time t = 0.2',
 '    Reading  Frame 5 at t = 0.25  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 5 at time t = 0.25',
 '    Reading  Frame 6 at t = 0.3  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 6 at time t = 0.3',
 '    Reading  Frame 7 at t = 0.35  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 7 at time t = 0.35',
 '    Reading  Frame 8 at t = 0.4  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 8 at time t = 0.4',
 '    Reading  Frame 9 at t = 0.45  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 9 at time t = 0.45',
 '    Reading  Frame 10 at t = 0.5  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 10 at time t = 0.5',
 '    Reading  Frame 11 at t = 0.55  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 11 at time t = 0.55',
 '    Reading  Frame 12 at t = 0.6  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 12 at time t = 0.6',
 '    Reading  Frame 13 at t = 0.65  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 13 at time t = 0.65',
 '    Reading  Frame 14 at t = 0.7  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 14 at time t = 0.7',
 '    Reading  Frame 15 at t = 0.75  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 15 at time t = 0.75',
 '    Reading  Frame 16 at t = 0.8  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 16 at time t = 0.8',
 '    Reading  Frame 17 at t = 0.85  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 17 at time t = 0.85',
 '    Reading  Frame 18 at t = 0.9  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 18 at time t = 0.9',
 '    Reading  Frame 19 at t = 0.95  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 19 at time t = 0.95',
 '    Reading  Frame 20 at t = 1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 20 at time t = 1.0',
 '    Reading  Frame 21 at t = 1.05  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 21 at time t = 1.05',
 '    Reading  Frame 22 at t = 1.1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 22 at time t = 1.1',
 '    Reading  Frame 23 at t = 1.15  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 23 at time t = 1.15',
 '    Reading  Frame 24 at t = 1.2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 24 at time t = 1.2',
 '    Reading  Frame 25 at t = 1.25  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 25 at time t = 1.25',
 '    Reading  Frame 26 at t = 1.3  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 26 at time t = 1.3',
 '    Reading  Frame 27 at t = 1.35  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 27 at time t = 1.35',
 '    Reading  Frame 28 at t = 1.4  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 28 at time t = 1.4',
 '    Reading  Frame 29 at t = 1.45  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 29 at time t = 1.45',
 '    Reading  Frame 30 at t = 1.5  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 30 at time t = 1.5',
 '    Reading  Frame 31 at t = 1.55  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 31 at time t = 1.55',
 '    Reading  Frame 32 at t = 1.6  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 32 at time t = 1.6',
 '    Reading  Frame 33 at t = 1.65  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 33 at time t = 1.65',
 '    Reading  Frame 34 at t = 1.7  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 34 at time t = 1.7',
 '    Reading  Frame 35 at t = 1.75  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 35 at time t = 1.75',
 '    Reading  Frame 36 at t = 1.8  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 36 at time t = 1.8',
 '    Reading  Frame 37 at t = 1.85  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 37 at time t = 1.85',
 '    Reading  Frame 38 at t = 1.9  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 38 at time t = 1.9',
 '    Reading  Frame 39 at t = 1.95  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 39 at time t = 1.95',
 '    Reading  Frame 40 at t = 2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 40 at time t = 2.0',
 '',
 '-----------------------------------',
 '',
 'Creating latex file...',
 "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ",
 '    already exists, files may be overwritten ',
 '',
 'Latex file created:  ',
 '  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/plots.tex',
 '',
 'Use pdflatex to create pdf file',
 'Created JSAnimation for figure 1',
 '',
 '--------------------------------------------------------',
 '',
 'Point your browser to:',
 '    file:///Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/_PlotIndex.html']
In [55]:
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[55]:


Once Loop Reflect

ODE Solver

In [116]:
rhol = 0.4
rhor = 0.9
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
In [52]:
%%system
mkdir _plots_ode_freeway2
Out[52]:
[]
In [117]:
Dt = tfinal/output_number
for i in range(output_number+1):
    clf()
    figure(figsize=(12,5))
    subplot(1,2,1)
    for j in range(0,N,4):
        plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
    plot(x[i][0:N:4],Dt*i*ones(N/4),'ob')
    plot([-N/rhol-10,tfinal+10],[Dt*i,Dt*i],'--r');
    axis([-N/rhol-10,tfinal+10,0,tfinal+1])
    title('Trajectory of Cars')
    xlabel('$x$')
    ylabel('$t$')
    
    subplot(1,2,2)
    plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
    axis([-N/rhol-10,tfinal+10,-0.1,1.1])
    title('Density')
    xlabel('$x$')
    ylabel('rho')
    fname = '_plots_ode_freeway2/frame%04i'%i + 'fig1' + '.png'
    savefig(fname)
<matplotlib.figure.Figure at 0x119bd9e50>
<matplotlib.figure.Figure at 0x119bd9410>
<matplotlib.figure.Figure at 0x10ffa7790>
<matplotlib.figure.Figure at 0x113f02110>
<matplotlib.figure.Figure at 0x1128c7e50>
<matplotlib.figure.Figure at 0x111674c50>
<matplotlib.figure.Figure at 0x10b9902d0>
<matplotlib.figure.Figure at 0x11a6c8150>
<matplotlib.figure.Figure at 0x11c6bc790>
<matplotlib.figure.Figure at 0x115305a90>
<matplotlib.figure.Figure at 0x10d6b9110>
<matplotlib.figure.Figure at 0x112955e50>
<matplotlib.figure.Figure at 0x113f47e10>
<matplotlib.figure.Figure at 0x1185644d0>
<matplotlib.figure.Figure at 0x10e133cd0>
<matplotlib.figure.Figure at 0x1116b4a10>
<matplotlib.figure.Figure at 0x1243cb110>
<matplotlib.figure.Figure at 0x113eecbd0>
<matplotlib.figure.Figure at 0x113f0b890>
<matplotlib.figure.Figure at 0x10ffbeb50>
<matplotlib.figure.Figure at 0x113f47d10>
<matplotlib.figure.Figure at 0x11ef0ac90>
<matplotlib.figure.Figure at 0x128d2e250>
<matplotlib.figure.Figure at 0x10e48d310>
<matplotlib.figure.Figure at 0x113ebb610>
<matplotlib.figure.Figure at 0x1128f3410>
<matplotlib.figure.Figure at 0x113f6bcd0>
<matplotlib.figure.Figure at 0x127572210>
<matplotlib.figure.Figure at 0x11166f2d0>
<matplotlib.figure.Figure at 0x11eefaa90>
<matplotlib.figure.Figure at 0x111772d50>
<matplotlib.figure.Figure at 0x10ec40d50>
<matplotlib.figure.Figure at 0x113e53050>
<matplotlib.figure.Figure at 0x11a73bdd0>
<matplotlib.figure.Figure at 0x10f739890>
<matplotlib.figure.Figure at 0x11ba2cd90>
<matplotlib.figure.Figure at 0x10f71bd10>
<matplotlib.figure.Figure at 0x10ec47d10>
<matplotlib.figure.Figure at 0x111650590>
<matplotlib.figure.Figure at 0x11a729e90>
<matplotlib.figure.Figure at 0x113f6a910>
<matplotlib.figure.Figure at 0x111674d10>
<matplotlib.figure.Figure at 0x11eefa810>
<matplotlib.figure.Figure at 0x1243c7490>
<matplotlib.figure.Figure at 0x11a767cd0>
<matplotlib.figure.Figure at 0x11ba0e450>
<matplotlib.figure.Figure at 0x111650ed0>
<matplotlib.figure.Figure at 0x10ecb8b90>
<matplotlib.figure.Figure at 0x10e107e50>
<matplotlib.figure.Figure at 0x10e12f710>
<matplotlib.figure.Figure at 0x10e48d410>
<matplotlib.figure.Figure at 0x11a73b210>
<matplotlib.figure.Figure at 0x11a7295d0>
<matplotlib.figure.Figure at 0x11baad2d0>
<matplotlib.figure.Figure at 0x1160ddc90>
<matplotlib.figure.Figure at 0x1185b6310>
<matplotlib.figure.Figure at 0x111763790>
<matplotlib.figure.Figure at 0x113e53710>
<matplotlib.figure.Figure at 0x10dd9f490>
<matplotlib.figure.Figure at 0x10e073450>
<matplotlib.figure.Figure at 0x1116092d0>
<matplotlib.figure.Figure at 0x11a6cf890>
<matplotlib.figure.Figure at 0x1153163d0>
<matplotlib.figure.Figure at 0x11e23e290>
<matplotlib.figure.Figure at 0x115305dd0>
<matplotlib.figure.Figure at 0x124386c10>
<matplotlib.figure.Figure at 0x11c61f190>
<matplotlib.figure.Figure at 0x128d1d750>
<matplotlib.figure.Figure at 0x1116cc3d0>
<matplotlib.figure.Figure at 0x11c6a8510>
<matplotlib.figure.Figure at 0x10ffbead0>
<matplotlib.figure.Figure at 0x120b93650>
<matplotlib.figure.Figure at 0x1115d5e10>
<matplotlib.figure.Figure at 0x10f65b1d0>
<matplotlib.figure.Figure at 0x11ba2c210>
<matplotlib.figure.Figure at 0x11c37f1d0>
<matplotlib.figure.Figure at 0x11e268410>
<matplotlib.figure.Figure at 0x10e4baa10>
<matplotlib.figure.Figure at 0x10dd90690>
<matplotlib.figure.Figure at 0x1140eb310>
<matplotlib.figure.Figure at 0x1160f0690>
In [118]:
figno = 1
fname = '_plots_ode_freeway2/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(12, 8),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[118]:


Once Loop Reflect

Case 3: $\rho_l<\rho_m<\rho_r$ and $\gamma/(1+\gamma)\geq \rho_l$

We construct the same convex hull as Case 2. There is only a shock in this solution.

$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_r \ \ \ \text{if}\ \ x > st, \end{array} \right.$$

where $s=\frac{f_\epsilon(\rho_r)-f_\epsilon(\rho_l)}{\rho_r-\rho_l}.$

In the numerical simulation, the initial condition we are using is:

$$\rho_l=0.2, \rho_r=0.9.$$

PDE Solver

In [59]:
%%system

make data ql=0.2 qr=0.9
make .plots
Out[59]:
['rm -f .data',
 'Wrote data to initial_data.txt',
 'python setrun.py              classic                  ',
 'touch .data',
 '/Library/Developer/CommandLineTools/usr/bin/make output',
 'rm -f .output',
 'python /Users/Peng/Documents/AMATH574/clawpack/clawutil/src/python/clawutil/runclaw.py xclaw                          _output                     \\',
 '\tTrue                    False                     . False False None',
 'Reading data file: claw.data   ',
 '         first 5 lines are comments and will be skipped',
 ' running...',
 '  ',
 'Reading data file: setprob.data',
 '         first 5 lines are comments and will be skipped',
 'CLAW1EZ: Frame    1 output files done at time t =  0.5000D-01',
 '',
 'CLAW1EZ: Frame    2 output files done at time t =  0.1000D+00',
 '',
 'CLAW1EZ: Frame    3 output files done at time t =  0.1500D+00',
 '',
 'CLAW1EZ: Frame    4 output files done at time t =  0.2000D+00',
 '',
 'CLAW1EZ: Frame    5 output files done at time t =  0.2500D+00',
 '',
 'CLAW1EZ: Frame    6 output files done at time t =  0.3000D+00',
 '',
 'CLAW1EZ: Frame    7 output files done at time t =  0.3500D+00',
 '',
 'CLAW1EZ: Frame    8 output files done at time t =  0.4000D+00',
 '',
 'CLAW1EZ: Frame    9 output files done at time t =  0.4500D+00',
 '',
 'CLAW1EZ: Frame   10 output files done at time t =  0.5000D+00',
 '',
 'CLAW1EZ: Frame   11 output files done at time t =  0.5500D+00',
 '',
 'CLAW1EZ: Frame   12 output files done at time t =  0.6000D+00',
 '',
 'CLAW1EZ: Frame   13 output files done at time t =  0.6500D+00',
 '',
 'CLAW1EZ: Frame   14 output files done at time t =  0.7000D+00',
 '',
 'CLAW1EZ: Frame   15 output files done at time t =  0.7500D+00',
 '',
 'CLAW1EZ: Frame   16 output files done at time t =  0.8000D+00',
 '',
 'CLAW1EZ: Frame   17 output files done at time t =  0.8500D+00',
 '',
 'CLAW1EZ: Frame   18 output files done at time t =  0.9000D+00',
 '',
 'CLAW1EZ: Frame   19 output files done at time t =  0.9500D+00',
 '',
 'CLAW1EZ: Frame   20 output files done at time t =  0.1000D+01',
 '',
 'CLAW1EZ: Frame   21 output files done at time t =  0.1050D+01',
 '',
 'CLAW1EZ: Frame   22 output files done at time t =  0.1100D+01',
 '',
 'CLAW1EZ: Frame   23 output files done at time t =  0.1150D+01',
 '',
 'CLAW1EZ: Frame   24 output files done at time t =  0.1200D+01',
 '',
 'CLAW1EZ: Frame   25 output files done at time t =  0.1250D+01',
 '',
 'CLAW1EZ: Frame   26 output files done at time t =  0.1300D+01',
 '',
 'CLAW1EZ: Frame   27 output files done at time t =  0.1350D+01',
 '',
 'CLAW1EZ: Frame   28 output files done at time t =  0.1400D+01',
 '',
 'CLAW1EZ: Frame   29 output files done at time t =  0.1450D+01',
 '',
 'CLAW1EZ: Frame   30 output files done at time t =  0.1500D+01',
 '',
 'CLAW1EZ: Frame   31 output files done at time t =  0.1550D+01',
 '',
 'CLAW1EZ: Frame   32 output files done at time t =  0.1600D+01',
 '',
 'CLAW1EZ: Frame   33 output files done at time t =  0.1650D+01',
 '',
 'CLAW1EZ: Frame   34 output files done at time t =  0.1700D+01',
 '',
 'CLAW1EZ: Frame   35 output files done at time t =  0.1750D+01',
 '',
 'CLAW1EZ: Frame   36 output files done at time t =  0.1800D+01',
 '',
 'CLAW1EZ: Frame   37 output files done at time t =  0.1850D+01',
 '',
 'CLAW1EZ: Frame   38 output files done at time t =  0.1900D+01',
 '',
 'CLAW1EZ: Frame   39 output files done at time t =  0.1950D+01',
 '',
 'CLAW1EZ: Frame   40 output files done at time t =  0.2000D+01',
 '',
 'Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL',
 '==> runclaw: Will take data from  /Users/Peng/Documents/AMATH574/am574-group07/presentation',
 '==> runclaw: Will write output to  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '==> runclaw: Removing all old fort files in  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '',
 '==> Running with command:',
 '    /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw',
 '',
 '==> runclaw: Finished executing',
 '',
 '==> runclaw: Done executing /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw via clawutil.runclaw.py',
 '==> runclaw: Output is in  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 '/Library/Developer/CommandLineTools/usr/bin/make plots',
 'rm -f .plots',
 'python /Users/Peng/Documents/AMATH574/clawpack/visclaw/src/python/visclaw/plotclaw.py _output                     _plots                     setplot.py           ',
 'Importing setplot.setplot from /Users/Peng/Documents/AMATH574/am574-group07/presentation.',
 'Executed setplot successfully',
 'Will plot 41 frames numbered: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]',
 'Will make 1 figure(s) for each frame, numbered:  [1]',
 '',
 '-----------------------------------',
 '',
 '',
 'Creating html pages for figures...',
 '',
 "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ",
 '    already exists, files may be overwritten ',
 'Now making png files for all figures...',
 '    Reading  Frame 0 at t = 0  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 0 at time t = 0.0',
 '    Reading  Frame 1 at t = 0.05  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 1 at time t = 0.05',
 '    Reading  Frame 2 at t = 0.1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 2 at time t = 0.1',
 '    Reading  Frame 3 at t = 0.15  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 3 at time t = 0.15',
 '    Reading  Frame 4 at t = 0.2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 4 at time t = 0.2',
 '    Reading  Frame 5 at t = 0.25  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 5 at time t = 0.25',
 '    Reading  Frame 6 at t = 0.3  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 6 at time t = 0.3',
 '    Reading  Frame 7 at t = 0.35  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 7 at time t = 0.35',
 '    Reading  Frame 8 at t = 0.4  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 8 at time t = 0.4',
 '    Reading  Frame 9 at t = 0.45  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 9 at time t = 0.45',
 '    Reading  Frame 10 at t = 0.5  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 10 at time t = 0.5',
 '    Reading  Frame 11 at t = 0.55  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 11 at time t = 0.55',
 '    Reading  Frame 12 at t = 0.6  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 12 at time t = 0.6',
 '    Reading  Frame 13 at t = 0.65  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 13 at time t = 0.65',
 '    Reading  Frame 14 at t = 0.7  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 14 at time t = 0.7',
 '    Reading  Frame 15 at t = 0.75  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 15 at time t = 0.75',
 '    Reading  Frame 16 at t = 0.8  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 16 at time t = 0.8',
 '    Reading  Frame 17 at t = 0.85  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 17 at time t = 0.85',
 '    Reading  Frame 18 at t = 0.9  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 18 at time t = 0.9',
 '    Reading  Frame 19 at t = 0.95  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 19 at time t = 0.95',
 '    Reading  Frame 20 at t = 1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 20 at time t = 1.0',
 '    Reading  Frame 21 at t = 1.05  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 21 at time t = 1.05',
 '    Reading  Frame 22 at t = 1.1  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 22 at time t = 1.1',
 '    Reading  Frame 23 at t = 1.15  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 23 at time t = 1.15',
 '    Reading  Frame 24 at t = 1.2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 24 at time t = 1.2',
 '    Reading  Frame 25 at t = 1.25  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 25 at time t = 1.25',
 '    Reading  Frame 26 at t = 1.3  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 26 at time t = 1.3',
 '    Reading  Frame 27 at t = 1.35  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 27 at time t = 1.35',
 '    Reading  Frame 28 at t = 1.4  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 28 at time t = 1.4',
 '    Reading  Frame 29 at t = 1.45  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 29 at time t = 1.45',
 '    Reading  Frame 30 at t = 1.5  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 30 at time t = 1.5',
 '    Reading  Frame 31 at t = 1.55  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 31 at time t = 1.55',
 '    Reading  Frame 32 at t = 1.6  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 32 at time t = 1.6',
 '    Reading  Frame 33 at t = 1.65  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 33 at time t = 1.65',
 '    Reading  Frame 34 at t = 1.7  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 34 at time t = 1.7',
 '    Reading  Frame 35 at t = 1.75  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 35 at time t = 1.75',
 '    Reading  Frame 36 at t = 1.8  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 36 at time t = 1.8',
 '    Reading  Frame 37 at t = 1.85  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 37 at time t = 1.85',
 '    Reading  Frame 38 at t = 1.9  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 38 at time t = 1.9',
 '    Reading  Frame 39 at t = 1.95  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 39 at time t = 1.95',
 '    Reading  Frame 40 at t = 2  from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output',
 'Frame 40 at time t = 2.0',
 '',
 '-----------------------------------',
 '',
 'Creating latex file...',
 "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ",
 '    already exists, files may be overwritten ',
 '',
 'Latex file created:  ',
 '  /Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/plots.tex',
 '',
 'Use pdflatex to create pdf file',
 'Created JSAnimation for figure 1',
 '',
 '--------------------------------------------------------',
 '',
 'Point your browser to:',
 '    file:///Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/_PlotIndex.html']
In [60]:
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[60]:


Once Loop Reflect

ODE Solver

In [119]:
rhol = 0.2
rhor = 0.9
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
In [62]:
%%system
mkdir _plots_ode_freeway3
Out[62]:
[]
In [120]:
Dt = tfinal/output_number
for i in range(output_number+1):
    clf()
    figure(figsize=(12,5))
    subplot(1,2,1)
    for j in range(0,N,4):
        plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
    plot(x[i][0:N:4],Dt*i*ones(N/4),'ob')
    plot([-N/rhol-10,tfinal+10],[Dt*i,Dt*i],'--r');
    axis([-N/rhol-10,tfinal+10,0,tfinal+1])
    title('Trajectory of Cars')
    xlabel('$x$')
    ylabel('$t$')
    
    subplot(1,2,2)
    plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
    axis([-N/rhol-10,tfinal+10,-0.1,1.1])
    title('Density')
    xlabel('$x$')
    ylabel('rho')
    fname = '_plots_ode_freeway3/frame%04i'%i + 'fig1' + '.png'
    savefig(fname)
<matplotlib.figure.Figure at 0x11a6c7750>
<matplotlib.figure.Figure at 0x11a6c7e10>
<matplotlib.figure.Figure at 0x10f62f250>
<matplotlib.figure.Figure at 0x11e21c590>
<matplotlib.figure.Figure at 0x11169bd90>
<matplotlib.figure.Figure at 0x1185b6690>
<matplotlib.figure.Figure at 0x118600cd0>
<matplotlib.figure.Figure at 0x120b1f6d0>
<matplotlib.figure.Figure at 0x1128f5890>
<matplotlib.figure.Figure at 0x11a6cffd0>
<matplotlib.figure.Figure at 0x11e23e890>
<matplotlib.figure.Figure at 0x11efc4610>
<matplotlib.figure.Figure at 0x113eec350>
<matplotlib.figure.Figure at 0x111797d10>
<matplotlib.figure.Figure at 0x11a6c71d0>
<matplotlib.figure.Figure at 0x10ff8dad0>
<matplotlib.figure.Figure at 0x113e53710>
<matplotlib.figure.Figure at 0x10e13b9d0>
<matplotlib.figure.Figure at 0x109f22dd0>
<matplotlib.figure.Figure at 0x10ec90510>
<matplotlib.figure.Figure at 0x11ba29090>
<matplotlib.figure.Figure at 0x11168ec10>
<matplotlib.figure.Figure at 0x114a5acd0>
<matplotlib.figure.Figure at 0x116114f10>
<matplotlib.figure.Figure at 0x11179ea10>
<matplotlib.figure.Figure at 0x1117474d0>
<matplotlib.figure.Figure at 0x1243cbf90>
<matplotlib.figure.Figure at 0x1117e4a50>
<matplotlib.figure.Figure at 0x10e0cab50>
<matplotlib.figure.Figure at 0x1129600d0>
<matplotlib.figure.Figure at 0x10f726d50>
<matplotlib.figure.Figure at 0x1117aba10>
<matplotlib.figure.Figure at 0x11a733650>
<matplotlib.figure.Figure at 0x124386c10>
<matplotlib.figure.Figure at 0x11a733ed0>
<matplotlib.figure.Figure at 0x1185f7bd0>
<matplotlib.figure.Figure at 0x113f0d750>
<matplotlib.figure.Figure at 0x113e65e10>
<matplotlib.figure.Figure at 0x1243b0d10>
<matplotlib.figure.Figure at 0x10d71bd90>
<matplotlib.figure.Figure at 0x11295d450>
<matplotlib.figure.Figure at 0x1116ae390>
<matplotlib.figure.Figure at 0x11166f390>
<matplotlib.figure.Figure at 0x11eff3290>
<matplotlib.figure.Figure at 0x1128ab390>
<matplotlib.figure.Figure at 0x10e107310>
<matplotlib.figure.Figure at 0x111654b50>
<matplotlib.figure.Figure at 0x11a733950>
<matplotlib.figure.Figure at 0x1128aefd0>
<matplotlib.figure.Figure at 0x11293a5d0>
<matplotlib.figure.Figure at 0x111650790>
<matplotlib.figure.Figure at 0x1116ae890>
<matplotlib.figure.Figure at 0x11a7108d0>
<matplotlib.figure.Figure at 0x1036b6f10>
<matplotlib.figure.Figure at 0x1116ae790>
<matplotlib.figure.Figure at 0x10eca8550>
<matplotlib.figure.Figure at 0x11a700b90>
<matplotlib.figure.Figure at 0x120b1e8d0>
<matplotlib.figure.Figure at 0x11f6244d0>
<matplotlib.figure.Figure at 0x1116a5e90>
<matplotlib.figure.Figure at 0x10e4eb890>
<matplotlib.figure.Figure at 0x124d84ed0>
<matplotlib.figure.Figure at 0x11e23ef10>
<matplotlib.figure.Figure at 0x1116af290>
<matplotlib.figure.Figure at 0x10e04e210>
<matplotlib.figure.Figure at 0x11a729ed0>
<matplotlib.figure.Figure at 0x11175e150>
<matplotlib.figure.Figure at 0x120b8c7d0>
<matplotlib.figure.Figure at 0x10f5f4c90>
<matplotlib.figure.Figure at 0x11169b350>
<matplotlib.figure.Figure at 0x1117ae590>
<matplotlib.figure.Figure at 0x1160dd090>
<matplotlib.figure.Figure at 0x111681850>
<matplotlib.figure.Figure at 0x120473ed0>
<matplotlib.figure.Figure at 0x11bab6390>
<matplotlib.figure.Figure at 0x128d2e3d0>
<matplotlib.figure.Figure at 0x11efcd890>
<matplotlib.figure.Figure at 0x1185a5b10>
<matplotlib.figure.Figure at 0x113e4fdd0>
<matplotlib.figure.Figure at 0x119bfe950>
<matplotlib.figure.Figure at 0x1128d8150>
In [121]:
figno = 1
fname = '_plots_ode_freeway3/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))

fig = plt.figure(figsize=(12, 8),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))


animation.FuncAnimation(fig, animate, init_func=init,
                              frames=len(filenames), interval=200, blit=True)
Out[121]:


Once Loop Reflect

Clustering of Night-time Driving Model

Introduction

When driving at night on an unfamiliar mountain road, facing with the empty road will cause the speed limitation by the uncertainty of the road condition. However if there are other cars in front of the road, it will be much easier to drive faster, since their tail lights indicate how the road twists and turns.

The velocity function for night-time driving is

$$ v(\rho)=\left\{ \begin{array}{l l} v_0 \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$

Hence,

$$ f(\rho)=\left\{ \begin{array}{l l} v_0\rho \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho^2 \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1\rho(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$

Clustering Behavior

We take uniformly spaced cars as initial density on an otherwise empty road.

In LeVeque's paper:

"The lead driver sees $\rho=0$ and drives at speed $v_0$, and then the second driver initially sees $\rho=\rho_0$ and can drive faster. But as the second car accelerates, the third driver observes a drop in density and hence this driver starts to slow down. This causes the fourth driver to go faster."

It has been proven that the uniformly spaced traffic with distance $\frac{1}{\rho_0}$ with speed $v(\rho_0)$ is unstable under pertubations.

Without Perturbation

In the car-following model, the parameters are: $$ \rho_a=0.1,\rho_b=0.3, c=10, v_0=1. $$

In [65]:
def vel_night(rho):
    rhoa, rhob = 0.1, 0.3
    v0, vmax = 1., 3.
    v1 = vmax/(1.-rhob)
    c = (vmax-v0)/(rhob-rhoa)
    n = rho.size
    v = zeros(n)
    for i in range(n):
        if rho[i] >= 0. and rho[i] < rhoa:
            v[i] = v0
        elif rho[i] >=rhoa and rho[i] <= rhob:
            v[i] = c*rho[i]
        elif rho[i] > rhob and rho[i] <= 1.:
            v[i] = v1*(1.-rho[i])
        else:
            print "Invalid density value!"
    return v
In [66]:
def ode_solver_night(N,tfinal,output_number,rhol,rhor):
    from copy import copy
    Dt = tfinal/output_number
    dt = min(0.01,Dt)
    x = zeros((output_number+1,N))
    rho = copy(x)
    x[0] = initialization(N,rhol,rhor)
    rho[0][1:N+1] = -1./diff(x[0])
    rho[0][0] = rhor
    m = int(Dt/dt)
    xs ,rhos = copy(x[0]), copy(rho[0])
    for i in range(1,output_number+1):
        for j in range(m):
            v = vel_night(rhos)
            xs = xs + v*dt
            rhos[1:N+1] = -1./diff(xs)
            rhos[0] = rhor
        v = vel_night(rhos)
        xs = xs + v*(Dt-dt*m)
        rhos[1:N+1] = -1./diff(xs)
        rhos[0] = rhor
        x[i], rho[i] = copy(xs), copy(rhos)
    return x, rho
In [123]:
N = 48
tfinal = 80
output_number = 80

figure(figsize=(15,5))
rhol = 1./5
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')

rhol = 1./4
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/4$')

rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
Out[123]:
<matplotlib.text.Text at 0x1160fe690>

With Perturbation

With random perturbations to the speeds, 2-car platoons can no longer be observed. Instead, there are random n-car platoons, which is more reasonable to model the night-time traffic.

In [69]:
def ode_solver_night_rand(N,tfinal,output_number,rhol,rhor):
    from copy import copy
    from numpy.random import rand
    A = 0.1
    Dt = tfinal/output_number
    dt = 0.01
    x = zeros((output_number+1,N))
    rho = copy(x)
    x[0] = initialization(N,rhol,rhor)
    rho[0][1:N+1] = -1./diff(x[0])
    rho[0][0] = rhor
    m = int(Dt/dt)
    xs ,rhos = copy(x[0]), copy(rho[0])
    for i in range(1,output_number+1):
        for j in range(m):
            v = vel_night(rhos) + A*(rand(N)-0.5*ones(N))
            xs = xs + v*dt
            rhos[1:N+1] = -1./diff(xs)
            rhos[0] = rhor
        v = vel_night(rhos) + A*(rand(N)-0.5*ones(N))
        xs = xs + v*(Dt-dt*m)
        rhos[1:N+1] = -1./diff(xs)
        rhos[0] = rhor
        x[i], rho[i] = copy(xs), copy(rhos)
    return x, rho
In [122]:
N = 48
tfinal = 80
output_number = 80

figure(figsize=(15,5))
rhol = 1./5
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')

rhol = 1./4
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/4$')

rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
    plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
Out[122]:
<matplotlib.text.Text at 0x1160e1250>

Reference

  • LeVeque, R. J. (2001). Some traffic flow models illustrating interesting hyperbolic behavior. Minisymposium on traffic flow, Chicago.
  • Wiens, J. K., Stockie, J. M., Williams, J. F. (2013). Riemann solver for a kinematic wave traffic model with discontinuous flux. Journal of Computational Physics, 242, 1-23.
  • Jiang, R., Wu, Q. S. (2007). The night driving behavior in a car-following model. Physica A: Statistical Mechanics and its Applications, 375(1), 297-306.
  • Dias, J. P., & Figueira, M. (2005). On the approximation of the solutions of the Riemann problem for a discontinuous conservation law. Bulletin of the Brazilian Mathematical Society, 36(1), 115-125.
  • Zhang, H. M. (2003). Anisotropic property revisited––does it hold in multi-lane traffic. Transportation Research Part B: Methodological, 37(6), 561-577.
In [ ]: